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Current renewed interest in exploration of the moon, Mars, and other planetary objects 
is driving technology development in many fields of space system design. In particular, 
there is a desire to land both robotic and human missions on the moon and elsewhere. The 
landing guidance system must be able to deliver the vehicle to a desired soft landing while 
meeting several constraints necessary for the safety of the vehicle. Due to performance 
limitations of current launch vehicles, it is desired to minimize the amount of fuel used. In 
addition, the landing site may change in real-time in order to avoid previously undetected 
hazards which become apparent during the landing maneuver. 

This complicated maneuver can be broken into simpler subproblems that bound the full 
problem. One such subproblem is to find a minimum-fuel landing solution that meets con- 
straints on the initial state, final state, and bounded thrust acceleration magnitude. With 
the assumptions of constant gravity and negligible atmosphere, the form of the optimal 
steering law is known, and the equations of motion can be integrated analytically, resulting 
in a system of five equations in five unknowns. It is shown that this system of equations can 
be reduced analytically to two equations in two unknowns. With an additional assumption 
of constant thrust acceleration magnitude, this system can be reduced further to one equa- 
tion in one unknown. It is shown that these unknowns can be bounded analytically. An 
algorithm is developed to quickly and reliably solve the resulting one-dimensional bounded 
search, and it is used as a real-time guidance applied to a lunar landing test case. 

I. Introduction 

Today is an exciting time for the exploration of the moon, Mars, and other planetary objects. Renewed 
interest in lunar exploration is driving new technology development for both human and robotic missions. 
NASA has begun a new era of space exploration with the Constellation Program, with a goal of returning to 
the moon to build a sustainable long term human presence. 1 The Lunar X-Prize has promised to reward $30 
million to the first private team to land a rover on the moon. 2 China, India, and Japan have all successfully 
orbited probes around the moon with hopes of someday putting down a lander. 3,4 While both Russia and 
the United States have previously landed vehicles on the moon and Mars, new technology development is 
necessary to achieve greater levels of performance and safety to successfully accomplish these missions. 

In the lunar landing scenario, the vehicle begins in some initial hyperbolic transfer orbit or a parking orbit 
around the Moon. Given the initial orbit and a desired landing point on the surface, a targeting algorithm 
determines the best Time-of-Ignition (TIG) at which to begin the deorbit. A burn is completed to put the 
lander in a coasting elliptical transfer orbit. At the designated time, the powered terminal descent phase 
begins. This phase can be broken into a braking phase and an approach phase. 5 During the braking phase, 
much of the velocity relative to the planet is reduced. During the approach phase, a pitch up maneuver 
places the vehicle in a vertical attitude above the landing site. The vehicle then descends vertically to the 
surface. During any part of the braking or approach phases, hazard avoidance may necessitate a target 
redesignation. 
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II. Current and Suggested Guidance Methods for Powered Landing 


II. A. Apollo Lunar Descent Guidance 

The heart of the Apollo lunar descent guidance is the assumption that the thrust acceleration is a quadratic 
function of time. This defines a two-point boundary value problem with five degrees of freedom which allows 
for the specification of the initial position, initial velocity, final position, final velocity, and final acceleration. 
These equations can be solved for the coefficients as a function of the known initial position and velocity, the 
desired final position, velocity, and acceleration, and the time-to-go ( t go )• The result is a simple guidance 
algorithm that produces feasible trajectories and which could run on the Apollo flight computer. In terms 
of fuel used, it has been shown that the Apollo guidance performs well when compared to a fuel optimal 
trajectory when only small target redesignations are necessary. However, when large target redesignations 
are needed, the Apollo guidance uses much more fuel than a fuel optimal trajectory. 6,7 Also, care must 
be taken to define a proper value of t go in order that the trajectory meet the desired characteristics. More 
details of the Apollo lunar descent guidance can be found in references 8, 9, and 10. 

II. B. Powered Explicit Guidance (PEG) 

The Powered Explicit Guidance (PEG) algorithm is used by the Space Shuttle for near fuel-optimal exoatmo- 
spheric powered maneuvers. It was derived for powered ascent trajectories which end in an orbital injection. 
As such, there is no constraint on the downrange position. Beginning with the optimal control law for pow- 
ered flight over a flat Earth in a uniform gravitational field, it is assumed that the unit thrust direction is 
mainly in the downrange direction (an assumption that holds for most Space Shuttle maneuvers). A further 
small angle assumption is made such that the unit thrust direction vector is approximated to be a linear 
function of time. With these assumptions, the equations of motion can be integrated analytically to produce 
a system of nine equations in nine unknowns. All variables can be computed given a value of the velocity-to- 
be-gained vector (V go ). A predictor-corrector algorithm is used to find the solution, given an initial guess for 
V go and the desired final velocity and position vectors (except the downrange component of position). 11-13 
Delporte and Sauvient reported an extension of the explicit guidance concept with constraints on the guid- 
ance law, such as altitude at first stage engine cut-off, thermal flux at first stage separation, and the landing 
footprint of the first stage. 14 Fill expanded the PEG equations to include a constraint on downrange. 12 This 
allows for a solution with a fully constrained final position and velocity, such as that of a powered terminal 
descent. Unfortunately, this expanded capability was not implemented in the PEG flight software. 

As can be seen, there are many approximations made which reduce the optimality of PEG. There is 
also no way to completely constrain the final position with the original implementation; however, Fill’s 
improvements should allow for this. In addition, while PEG is able to combine multiple thrust arcs, such as 
a constant thrust arc followed by a constant acceleration arc, there is no computation of optimal switching 
times for the throttle command between maximum and minimum thrust levels. In fact, PEG was derived 
assuming that the thrust profile is a known function of time. 

II. C. Other Methods 

These methods have been suggested but have never flown on an actual space vehicle. 

II.C.l. Analytical Near Fuel- Optimal 

Since the fuel-optimal problem is difficult to solve, it has been proposed to solve related problems that are 
near fuel-optimal. D’Souza developed a guidance algorithm for powered vacuum flight in a uniform grav- 
itational field by minimizing a weighted function of the time-of-flight and the square of the acceleration 
magnitude. 15 Two similar ideas are proposed by Najson and Mease: 1) minimize the integral of the sum of 
the absolute values of the thrust acceleration components, and 2) minimize the square of the thrust accel- 
eration magnitude. 16 Application of optimal control theory results in an analytical solution in terms of the 
initial position/velocity state, desired final position/velocity state, and time-of-flight. D’Souza analytically 
calculated the optimal time-of-flight as a function of the state, while Najson and Mease used a numerical 
search to find the time-of-flight with the smallest fuel expenditure. For an engine with a constant specific im- 
pulse (I S p, Ross has shown that a fuel optimal solution is found by minimizing the integral of the magnitude 
of the thrust acceleration. 17 Thus, while these methods are fuel efficient, they are not fuel optimal. 
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Ueno and Yamaguchi derived an analytical guidance law by making several approximations . 18 Starting 
with the problem of fuel-optimal planar flight in a uniform gravitational field with a constant thrust accel- 
eration, a small angle assumption is made that assumes the flight is mostly horizontal with the thrust being 
directed nearly opposite the velocity vector. Constraints are placed on the final velocity vector and the final 
altitude, and the final crossrange and downrange are free. Optimal control theory then leads to a polynomial 
guidance law where the downrange and crossrange velocities are linear functions of time, the altitude rate is 
a quadratic function of time, and the altitude is a cubic function of time. This guidance algorithm is shown 
to work well for the defined problem. However, it does not allow for designation of a specific landing site, 
except by proper selection of the deorbit state. This fact and the assumption of near horizontal flight implies 
that the proposed algorithm will not work for hazard avoidance maneuvers. The approximations also imply 
that the problem is not truly fuel optimal. 

Uchiyama used a barrier function method to transform the fuel optimization problem with bounded thrust 
magnitude into an unconstrained problem . 19 The controls are taken as the total acceleration components 
along the radial and local horizontal directions, including both the thrust and gravitational acceleration 
terms. The magnitude of the total acceleration is then minimized. In order to find the necessary thrust 
acceleration commands, the gravitational vector is subtracted from the optimal total acceleration vector. 
Issues with this approach include the lack of a constraint on the thrust acceleration magnitude and a lack 
of fuel optimality. 

II. C. 2. Analytical Gravity Turn 

Mclnnes and Chornel developed guidance algorithms based on the gravity turn where the thrust vector 
is directed opposite the velocity vector . 20,21 The gravity turn has the desirable properties that the final 
velocity magnitude is zero and the final vehicle attitude is vertical. For planar motion and constant thrust 
acceleration, it is possible to analytically integrate the equations of motion for altitude, downrange, velocity 
magnitude, and time as a function of flight path angle. The trajectory can be shaped by breaking it into 
segments with different thrust acceleration magnitudes. This guidance algorithm is shown to work well, and 
Cliomel has demonstrated that it can be used for target redesignations. However, while the gravity turn may 
be near fuel-optimal for descent from orbit, it is not a fuel optimal maneuver for target redesignations near 
the end of the terminal descent trajectory. Thus, w r hile this method is a viable solution for an analytical 
terminal descent guidance, it may not be applicable to all terminal descent situations. 

II. C. 3. Numerical Methods 

Acikmese and Ploen have suggested using a direct numerical method to solve the powered terminal descent 
problem . 22 The problem of time-fixed, fuel-optimal flight in a uniform gravitational field is made convex 
under certain assumptions. This convex problem is then converted into a second-order cone programming 
problem (SOCP) and solved using a numerical interior point method to find a global optimum with a known 
upper bound on the number of iterations. An outer iteration is necessary to find the optimal time-of- flight. 
The proposed guidance algorithm would numerically solve the SOCP onboard the vehicle. This method has 
the advantage of allowing several constraints, such as bounded thrust, bounded vehicle attitude during flight, 
constrained final vehicle attitude, and constraints on altitude such that the vehicle does not fly subsurface. 
While this approach holds much potential, it has the disadvantage of using a numerical solver. Also, the 
outer iteration to find the best time-of-flight is undesirable since the SOCP must be solved multiple times. 
While a solution to the SOCP can theoretically always be found, it would be difficult to verify convergence 
for all possible scenarios. Cognizant of these issues, a method was proposed whereby a family of optimal 
trajectories would be generated pre-flight using this method, and a table look-up would be used onboard to 
generate near fuel-optimal trajectories for any initial state . 23 

II. D. Areas for Improvement and Scope of Work 

A fuel-optimal guidance algorithm is needed for hazard avoidance and precision landing. It can be seen that 
many methods for powered terminal descent have been proposed. Only the Apollo Lunar Descent Guidance 
and the Shuttle Powered Explicit Guidance have been proven in flight. Most of the methods either make no 
use of optimal control theory, make approximations of the fuel optimal control law, or minimize performance 
functions that do not lead to fuel optimal solutions. Current mission designs put a high value on any design 
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feature that can reduce mass, thus a fuel optimal solution is desirable. Most of the methods also do not have 
the ability to bound the pitch angle of the vehicle during flight. This capability may be required to ensure 
that the vehicle attitude allows astronauts or sensors to view the landing site through a window. Some of 
the methods can be used for hazard avoidance maneuvers and some can not. 

It is desired to find a guidance algorithm which is fuel-optimal while meeting constraints on the final state, 
thrust magnitude, sensor and/or window viewing angles, and minimum approach altitude. The guidance 
must be able to generate commands in real-time without relying on the solution from the previous guidance 
cycle, as in the case of a target redesignation in mid-flight due to hazard avoidance at the original landing 
site. It is likely that this challenging problem has no analytical solution. However, it can be broken down 
into simpler subproblems which can be solved analytically. The solutions to these subproblems may be 
useful bounds for a numerical solution of the full problem. Towards that goal, this paper investigates the 
simpler subproblem of fuel optimal flight with bounded thrust acceleration magnitude from a given initial 
state to a desired final state; the pitch up maneuver to place the lander in the desired final vertical attitude 
is not considered. It will be shown that this problem can be reduced to an analytically bounded search with 
excellent convergence properties. 


III. Problem Definition 


The following assumptions are made: 

1. Atmospheric forces can be neglected. 

2. Gravitational acceleration is constant. 

3. The rotation of the planetary object can be neglected. 

4. The vehicle carries a perfectly expanded chemical rocket engine with which to create thrust force. 

5. The nozzle exit velocity of the propellant is a known constant. 

6. The thrust magnitude has known limits. 

7. The thrust direction is along the roll axis of the vehicle and can be commanded in any direction 
inst ant aneously. 

8. The final time is free. 


The only forces acting on the vehicle are the force due to gravity and its own thrust force. The equations 
of motion are: 


R = V 

V = a T u + g (1) 


where R is the position vector, V is the velocity vector, g is the gravitational acceleration vector, a-r is the 
magnitude of the thrust acceleration, and u is the unit thrust direction vector. 

It is desired to minimize the propellant used during the maneuver. It can be shown that minimizing 
the amount of propellant burned is equivalent to minimizing the integral of the magnitude of the thrust 
acceleration. 24 Note that no assumption of constant thrust or constant mass is made. 


{min J = m 0 — m /} ^ 


■| min 



( 2 ) 


where m 0 is the initial vehicle mass and rrif is the final vehicle mass. 

There is a maximum rate at which propellant mass can be expelled from the chemical rocket engine. The 
absolute minimum rate would be zero. However, there may be operating conditions such that the minimum 
rate is nonzero. These limits on propellant mass flow rate translate directly into limits on thrust magnitude. 
Note that it is not assumed that the thrust magnitude limits are constant. 


0 < T min (t) < T < T max (t) 


( 3 ) 
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Dividing equation 3 by the mass and substituting ot = yields 


0 < 


T min ( t ) 


< CLt ^ 


(t) 


m (t) in ( t ) 

The thrust acceleration magnitude limits are rewritten as: 

0 < aT min (t) < ar < ar, 

Constraints are placed on the final state such that: 


< 


(t) 


( 4 ) 


( 5 ) 
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= 0 
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= 0 

a Vy (Vy f 

-v Vft ) 

= 0 

a Vz (V Zf 


= 0 


( 6 ) 


where 

<Ji = on/off switch for each constraint (=0 -OR- 1) 

Rf = final position vector 

Vf = final velocity vector 

Rf g = specified final position vector 

Vf e = specified final velocity vector 

The <7i are added to allow easy comparison of solutions with different final constraints. 

IV. Application of Optimal Control Theory 

The problem can be solved with the well known techniques of optimal control theory, and it is a very 
well solved problem. As such, the results will be quoted directly here. Details of the solution can be found 
in references 24-33, as well as others. In this research, the terminology and methodology from reference 34 
are used. An assumption of familiarity with these methods is assumed. The Hamiltonian is given by: 


H = clt + A r ■ V + Ay • [cltu + g\ 


( 7 ) 


where Xp> and Ay are the costate Lagrange multipliers associated with position and velocity, respectively. 
Define the following variables: 


( 8 ) 



’ c Rx ' 
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Cn y 

= 

V Ry a Ry 


. C 'R . . 


VR Z VR Z 

II 
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Cy 

v y 

= 
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. C Vz _ 


V v z &v z 


(9) 


The Lagrange multipliers are continuous functions and can be written as: 24 

A r = Cr 
Ay = 


The total time-of-flight is defined as 


C v + C r t 

(10) 

as 


(t f - t ) 

(11) 

( tf - 1 0 ) 

(12) 
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The length of Ay is given by 

Ay = 

Ay 

= \JD + 2 Et + Ft 2 

(13) 

where 


D 

= CyCv 




E 

= Cr-C v 




F 

= Cr-Cr 

(14) 


IV. A. Bilinear Tangent Steering 


Lawden was the first to have derived the linear tangent steering law for powered vacuum flight in a uniform 
gravitational field. 25,26 Miele showed that the most general optimal control law for planar powered vacuum 
flight in a uniform gravitational field is the bilinear tangent steering law. 2 ' In vector notation, the optimal 
unit thrust vector in a uniform gravitational field is opposite the direction of the Ay Lagrange vector, or 
’’primer” vector: 24 



(C v + C r t ) 
sjD + 2 Et + Ft 2 


(15) 


If the motion is restricted to the XZ-plane and the angle of the unit thrust vector measured from the 
local horizontal, 0, is found, then equation 15 reduces to the more familiar bilinear tanget form: 


tan 9 = 


C Vz + Cr z t 

Cv x + Cr x t 


(16) 


Note that the structure of the bilinear tangent steering law is independent of the thrust function; the 
form of the bilinear tangent steering law holds for any bounded thrust profile. It is also independent of 
the optimization function. However, this steering law can reduce to various forms, such as linear tangent or 
constant thrust direction, depending on the optimization function, initial constraints, and final constraints. 27 
For example, when the downrange component of position is unconstrained, the bilinear tangent law reduces 
to a linear tangent law. 28,29 

The derivation of the classic bilinear tangent steering law uses the assumption of a uniform gravitational 
field. The gravitational acceleration can be reasonably assumed constant when the change in radial distance 
from the center of the planetary object during the maneuver is small. Yang found that when the bilinear 
tangent steering law is used in an inverse square gravitational field, it “is an excellent approximation to the 
true optimal control for powered flight arcs up to at least 15 degrees” of central arc. 35 In addition, it is 
possible that this approximation may hold reasonably well for trajectories with larger central arcs. 


IV. B. “Bang-Bang” Thrust Profile 

Leitmann showed that for a bounded thrust problem, the optimal thrust profile consists of up to three 
maximum or minimum thrust arcs. The structure of the thrust profile is such that the possible arcs include: 
max, max-min, max-min-max, min-max, or min. He showed that the use of intermediate thrust arcs is 
not optimal in most situations. In addition, he showed that the Hamiltonian, H , and costates, A r and 
Ay, are continuous across the thrust switches. 30,31 This means that the bilinear tangent steering law is a 
continuous function even when the thrust magnitude switches instantaneously. A solution to the special case 
of pure vertical landing was found by Meditch. He showed that the optimal thrust profile would be either 
max or min-max. 36 These thrust profiles are called “bang-bang” because the thrust magnitude “bangs” 
instantaneously between its maximum and minimum possible values. Note that modern rocket engines can 
be throttled very quickly. Topcu, Casoliva, and Mease noted that the engines for the Mars Smart Lander, 
which is to be launched in 2011, can be “throttled between minimum and maximum thrust in 30-40 ms, 
and so the thrust magnitude dynamics are very fast compared with the translational and mass dynamics” . 32 
This means that guidance commands that call for instantaneous thrust switches between maximum and 
minimum limits will not pose a problem for modern planetary landers. 
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The thrust acceleration magnitude switches between its upper and lower limits according to the value of 
the thrust acceleration magnitude switching structure (SS). 24 


S', S' = - 




\JD + 2 Et + Ft 2 


(17) 


Recall that the Lagrange multipliers Ay are continuous. Therefore, the switching structure is continuous. 
The thrust acceleration magnitude switching condition is given by 


SS < 1 

when 

CLT = &T m in 


SS > 1 

when 

CLt = a T mELX 


1— 1 
II 

co 

co 

when 

CLt < CLT < CLt 

-l min - 1 1 max 

(18) 


One question to ask is when an optimal, off-boundary acceleration subarc can occur. It is known that 
this can occur only when SS = 1. By examining the form of the SS, it can be seen that off-boundary 
acceleration subarcs will generally occur only instantaneously. However, there is a special case when the 
thrust acceleration can be continuously off-boundary for an entire trajectory. Examination of equation 17 
shows that a continuous off-boundary acceleration subarc can occur if 


D — Cy • Cy = 1 
E = Cr - Cy = 0 

F = C R -C R = 0 (19) 


This implies that a continuous off-boundary acceleration subarc can occur when Cr = A r = 0. Thus, for an 
off-acceleration boundary subarc: 

Ay = Cy = constant vector (20) 

Recall that the Lagrange multiplier Ay is continuous on an optimal trajectory. This implies that it has 
the same constant value for the entire trajectory. Thus, for an off-acceleration subarc to occur, the entire 
optimal trajectory must be off the acceleration boundaries. This means that the unit thrust direction vector 
will have the same constant value for the entire optimal trajectory. Thus, for an off-acceleration boundary 
subarc, we have 

u = constant unit vector (21) 


Recall that when SS=1, this means that 


Ay 


= y/D + 2 Et + Ft 2 = 


thrust acceleration is off-boundary for an entire subarc, Cr = A r = 0. 
Hamiltonian in equation 23 becomes 

H = Ay ■ g = Cy ■ g = 0 


1. For the special case when the 
Also, Ay = constant. Thus, the 

(22) 


This implies that Ay must be perpendicular to g. This implies that the unit thrust direction is in the 
horizontal plane. 

From equation 8, Cr is zero when ctr is zero. Recall from equation 6 that the a r are binary variables 
that act as on/off switches for the final position constraints. Thus, there are two cases when an off-boundary 
thrust acceleration solution is possible. The first case is when no final position is specified. The second case 
is when the desired final position is achieved by chance without a specific constraint on the final position. 
In this case, the thrust directions will only allow for a very limited set of trajectories. Thus, while this 
case should not be ignored, it will almost never appear in most real-world problems if the full final state is 
specified. Thus, most problems will have a thrust acceleration on one of its boundaries with instantaneous 
switches between the boundaries. However, it should be noted that although it is only possible to have 
an off-boundary acceleration solution when there is no final position constraint, it is also possible to have 
an on-boundary acceleration solution when there is no final position constraint. However, an on-boundary 
acceleration solution with no final position constraint will not have the condition of equation 22. This is 
discussed further in Section XI. 
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IV. C. Integrals of the Problem 

Certain first integrals of the optimization problem are known. 


• Hamiltonian 

For a time independent gravitational field, the Hamiltonian is constant. In addition, if the final time 
of the problem is free, the Hamiltonian has a constant value of zero. 33 Substiting for the optimal unit 
thrust direction into equation 7, the first integral of the problem is then 


H 

H 


CLt 

CLx 


1 - 


Ay 


+ Xr ■ V + Xy ■ g = 0 


1 - y/D + 2 Et + Ft 2 ] + (c R • v) + (c v + C r t ) • g = 0 


(23) 


• Rotational Symmetry 

The optimal powered flight problem has an axis of symmetry. Problems in an inverse gravitational 
field have a spherical symmetry. This leads to the following relationship: 

R x Xr + V x Xy = Constant (Inverse gravitational field) (24) 


A derivation of this relationship can be found in reference 33 (pp. 85-86). For problems with a uniform 
gravitational field, the axis of symmetry is about the gravitational vector. Thus, the component of the 
vector expressed in equation 24 that is parallel to the gravitational vector is constant. 32,33 


R X Xr + V X Ay 


• g = Constant 


(Uniform gravitational field) 


(25) 


V. Analytical Integration of the Equations of Motion 

It is possible to analytically integrate the equations of motion with the bilinear tangent steering law 
defined in equation 15 using standard integral tables. Yang integrated the equations with constant thrust 
magnitude. 3 ' Fill integrated them assuming that the thrust acceleration profile is a quadratic function of 
time. 12 Delporte and Sauvient integrated the equations for a constant thrust and the assumption that the 
two vectors defining the bilinear tangent law are perpendicular (Cr ■ Cy = 0). 14 Feeley integrated both 
the linear tangent guidance law constrained to a plane and the three dimensional bilinear tangent guidance 
law, both with constant thrust magnitude. 38,39 Leung, Calise, and Hull integrated the equations of motion 
constrained to a plane with the bilinear tangent law and constant thrust magnitude. 40,41 

For all these integrations, the result is a system of nonlinear equations for the final state as a function of 
the vectors defining the bilinear tangent law and the time-of- flight. This system of equations is difficult to 
solve, and in all the references where analytical integration was completed, the resulting system of nonlinear 
equations is solved using numerical methods. Indeed, Yang himself expressed that once the equations of 
motion are integrated analytically, “the determination of the constants involves the solving of 15 coupled 
nonlinear simultaneous equations. It is impractical to do it analytically. With the assistance of computers, 
a numerical scheme is much more appealing” , 37 However, it will be shown that analytical reduction of the 
problem is possible and yields useful results. 

It is convenient to integrate the equations of motion in terms of the time-to-go, r, defined in equation 
11. This change of variable is equivalent to integrating the equations of motion backwards from the final 
time. Consider the integration of the trajectory from the final time of flight to the initial time of flight. At 
t = tf, Tf = 0, V = Vf and R = Rf. At t = t Q , t = t q = ( tf — t a ), V — V 0 and R = R 0 . After completing 
the integrations and rearranging the equations, the following vectors can be defined: 



It is known that the optimal thrust acceleration magnitude will have up to three phases, with the thrust 
acceleration magnitude equal to either its minimum or maximum value. Only in certain very rare situations, 
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it is possible for ot to have a continuous, time-varying, off-boundary solution with a constant unit thrust 
direction. These cases will not be considered further. At this point, it will be assumed that the thrust 
acceleration limits are constant (i.e. (iT min =constant, &T max =constant). This means that ap will have a 
constant value on each phase. It can be shown that the same general solution exists for optimal trajectories 
with one, two, or three thrust acceleration subarcs, that is 


V 

R 


Pi 

F 


EG '\r + ( G ' 

-7F>) Cr + \7¥ 

, E 2 D 

! 7F 


Cy 


Vf*) f F 2 


<?*+(§- 


EGA 

Vf*J 


Cv 


(27) 


These equations represent the general integrated equations of motion for optimal solutions. They apply 
regardless of how many thrust acceleration subarcs exist (one, two, or three). Pi, Gi, and K, are integration 
constants where the subscript i refers to the number of subarcs. Details of the integration can be found in 
reference 24. 


VI. Thrust Vector Plane 

Lawden observed that the optimal thrust vector from equation 15 lies in an inertial plane defined by 
the two constant vectors Cr and Cy 26 (pp- 70-71). It never leaves this plane. From equation 27, it can 
be seen that the two vectors V and R can be written as linear combinations of Cr and Cy. Thus, all four 
vectors are coplanar. Define this plane as the Thrust Vector Plane, or simply the Thrust Plane (TP). A 
coordinate system can be set in this plane with the Z-axis perpendicular to the thrust plane. Marec observed 
that the vectors Cr and Cy can be written in terms of a magnitude and an angle 29 (p. 72). A seemingly 
new observation, but one that is obvious, is that the thrust plane is defined entirely by V and R, which are 
functions of the initial state, final state, gravity, and time-of- flight. This fact can be used to greatly simplify 
the problem. 

Define a cartesian coordinate system such that the Z-axis points up along the vertical direction with the 
origin at an altitude of zero. The X-axis and Y-axis form the horizontal plane. This coordinate system 
will be referred to as the inertial frame and will be denoted with a subscript of “I”. A subscript of “TP” 
will be used to denote vectors on the Thrust Plane. Define the Thrust Plane X-axis along the direction of 
Vj. Define the Thrust Plane Z-axis as a unit vector perpendicular to the plane. The Thrust Plane Y-axis 
completes the right-handed coordinate system. The axes of the Thrust Plane Coordinate System are then 
defined in terms of inertial vectors as 


XTP i 

VTPt 

ZtPj 


Vi 

Vi 

ZPP ! X XTPj 
V, x R, 

V/ x Ih 


(28) 


Define the angle 4> as the angle of the vector of interest off of the Thrust Plane. Define the angle 9 as 
the angle in the plane between the vector of interest and the V vector (i.e. Thrust Plane X-axis), measured 
positive about the Thrust Plane normal. The vectors R 0 , V 0 , Rf, Vf, and g become: 


R 0 

— Ro 

Vo 

= 

Rf 

= Rf 

Vf 

= v f 

9 

= 3 (' 


3 <t> Ro cos ®R 0 iTP + cos <j>R o sin 9 Ro j TP + sin (f>R 0 kpp 
! (f>y o cos 9y o ipp + cos (j)y o sin 9y o jp p + sin <j)y o kpp^j 
s <j>R f cos 0 Rf i T p + cos (f>R f sin 9R f jTP + sin (j> Rf kpp 
i 4>y. cos Qy. ipp + cos <j)y f sin 9y, jp p + sin ( 
g (cos (j> g cos Qgipp + cos <p g sin 9 g jpp + sin 4> g kpp 


Kpp 


(29) 
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Note that the vectors Cr, Cy, R, and V lie on the Thrust Plane. Thus <j> = 0 for each vector, and they 
have no Z-component. So, on the Thrust Plane we have 



V TP = V ma g (cOS Qyipp + Sin dyjpp'j 
Rtp = Rrnag (cOS 9r}PP + SUl OrJtP^J 



Cv T p = C y (cosOcyiTP + sm9c v jTP^ 



Cr tp = Cr (cos dc R irp + sin 0c R jpp^ 

(30) 

where 

Vmag = Magnitude of the V vector 
Rmag = Magnitude of the R vector 
G y = Magnitude of the Cy vector 

Cr = Magnitude of the Cr vector 


Note that 0y = 0 by definition, because the X-Thrust Plane axis is in the direction of V. However, for 
purposes of a more general derivation, this fact will be ignored for now. Define 


4> := (cos 0 Cr cos 9 Cv + sin 6 Cr sin 6 Cv ) 

(31) 

Note that this is 
equations 14 as 

equivalent to the cosine of the angle between Cy and Cr. 

Using these definitions, rewrite 


D = Cy • Cy = Cy 
E = Cr-Cv=CrCv^ 
F = Cr ■ Cr = C 2 r 

(32) 


VII. Integration Constants for Optimal Trajectories 

The integration constants (Pj, Gi , K, ) for each possible trajectory are defined using Thrust Plane nota- 
tion. 


VII. A. One-Phase Solution: ap = constant 

P, G, and I\ for a solution with one thrust acceleration subarc are given by: 24 


Pi = ar 


C v -\/c v 2 + 2 C r Cv^>t 0 + C r 2 t* 


C v 2 + 2 CrC v <S>t 0 + C r 2 t2 + C r t 0 + C y $ 


Gl = ' aTln ' C y + G y $ 

K\ = -O-tTo-J Cy 2 + 2 CrCv^t 0 + C_r 2 t 2 


(33) 


VII. B. Two-Phase Solution: O^T not = a Tmax 5 f i n = 

P, G, and K for a solution with two thrust acceleration subarcs with ap = clt , clt, = &p are given 
by: 24 
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= a Tmin (Cy - 1) - a Tmax (y/Cv 2 + 2 C r C v *t 0 + C r 2 t* - l) 


= -In 


- In 


1 + y] C v 2 $ 2 - ( Cy 2 - 1) 

Cy + Cy<I> 


!c v 2 + 2C fl Cy$r 0 + C fl 2 r 2 + C*r 0 + Cy4> 
1 + */Cy 2 $ 2 - (Cy 2 - 1) 


= 2 " a r m ax) 


-Cy$ + JCy 2 $ 2 - (Gy 2 - 1) 


+ 7j CL T Ia a. x T o\J Cy 2 + 2C R Cy$T 0 + Cfl 2 T 2 (34) 

VII. C. Two-Phase Solution: CL T no t = a Tmi n ? O'Tfin, = UTmax 

P, G, and K for a solution with two thrust acceleration subarcs with clt t = a-T ■> a T r = flr are given 

by .24 

p 2min - max = a T ma x ( C V ~ 1) “ ^T min (^J Cy 2 + 2 C R Cy$T 0 + C R T% - 1^ 


= — In 


— In 


1 - JCy 2 & - (Gy 2 - 1) 


Gy + Gy<l> 


! Cy 2 + 2GflGy$r 0 + Cr 2 t% + C R r 0 + Gy$ 
1 - ^Cy 2 * 2 - (Gy 2 - 1) 

/-Gy$ - Jc V 2 & - (Gy 2 - l) 


= 2 ( flT “ “ a T min ) 


+ 2 a T miTl T o\l Cy 2 +2C R Cy<&T 0 + Cr 2 T% 


VII. D. Three-Phase Solution 

P, G, and A' for a solution with three thrust acceleration subarcs are given by: 2 ^ 

P 3 = «T m „ Gy - Gy 2 + 2GflGy$r 0 + Cfl 2 T 2 

[ / v/ Gy 2 + 2G fl Gy$r 0 + C fl 2 r 2 + G fl r 0 + Cy4> 


G, = - In 


— In 


Gy + Gy$ 


1 - Y Gy 2< l> 2 - (Gy 2 - 1) 

1 + Jcv^^iC^ 2 ^) 


( flT max flT min ) 


K3 = 2 aT max r ° Y^ Cv 2 + 2 G fl Gy$r 0 + Gfl“r 2 

^Cy 2 ® 2 - ( Gy 2 - 1) 
- ( a r max - °T min ) - 77 
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VIII. Coefficient Relationships 


Starting with equations 27 and applying various combinations of vector dot products and vector cross 
products with Cr, Cy, R, and V, it is possible to derive the following four relationships. 


P 

G 

K 

0 


V 


( C R 

Vf [f (c v -V)-e(c r - h ) 

(. DF - E 2 ) 


-(S r .r)-\(c v .v) 

{Cr x Rj + (C'V x v'j 


(37) 


The derivations are straight-forward but cumbersome. For details, see Appendix C of reference 24. 

Equations 37 represent a system of four equations in five unknowns (Cy, Cr , 9c v , 9c R , and r Q ). The first 
of equations 37 is equivalent to the condition that the Hamiltonian is constant along the optimal trajectory 
from equation 23. It can be derived by subtracting the final value of the Hamiltonian from the initial value 
of the Hamiltonian. The fourth of equations 37 is equivalent to the condition of rotational symmetry from 
equation 25. It can be derived by subtracting the condition applied at the initial time from the condition 
applied to the final time. The fifth equation to complete the system comes from the knowledge that the 
Hamiltonian has a value of zero. It will be convenient to apply this condition at the final time. 


H f = a Tf 


i -Vd\ 


+ Cr ■ Vf + Cy ■ g = 0 


In Thrust Plane notation, equations 37 become 

P = C R V mag (cos 9 Cr cos dy + sin 9 Cr sin 9y ) 

r - —V ( cos sin 8y - siu Qc r cos 9y) 
Cy ma9 (cos 9c R sin 9c v - sin 9c R cos 9c v ) 
K = —C R R mag (cos 9 Cr cos 9 r + sin 9 Cr sin 9 R ) 

~ \ Cy Vmag (cos 9 Cv cos 9y + sin 9 Cv Sin 9y ) 

0 = CRRmag (cos 9 Cr sin 9ft - sin 9 Cr cos 9^) 

+ Cy V mag (cos 9 Cv sin 9y - sin 9 Cv cos 9y ) 


and the final Hamiltonian condition is 


(38) 

(39) 

(40) 

(41) 

(42) 


0 = a-jy [1 — Cy\ 

+ Cr Vf cos <j>v f (cos 9c R cos 9y f + sin 9c r sin 9y f ) 

+ Cy g cos (j) g (cos 9c v cos 9 g + sin 9c v sin 9 g ) = 0 (43) 


Equations 39 through 43 represent a system of five equations in five unknowns: Cy, Cr, 9c v , 9 c r , and r Q . 


IX. Reduction of Optimal Problem to Two Variables 


It is possible to analytically reduce the optimal problem to two variables. For these derivations, it is 
assumed that a constant thrust direction solution is not valid. This is equivalent to the assumption that all 
of the vectors Cr, Cy, R, and V are non-zero and that none of them are parallel. Recall that the Thrust 
Plane Coordinate system X-axis was defined as along the direction of V . Thus, the angle 9y is equal to zero 
by definition. However, for purposes of a general derivation, this fact will not be used. The first step is to 
solve equations 42 and 43 for Cr and Cy in terms of 9c R , 9c v , and r Q . Start with equation 42 and solve for 
Cr to obtain 


Cr 


— Cy ^ ma 9 


Rn 


(cos 9c v sin 9y 
(cos 9c r sin 9fi 


sin 9c v cos Gy) 
sin 9c r cos Off) 


(44) 
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Now take the final Hamiltonian in equation 43 and solve for Cr 


C r = - 


ctT f [1 — Cy\ + Cyg cos <j)g [cos 8c v cos 6 g + sin 0c v sin 0 g \ 


Vf cos cj>v f [cos 9c R cos 6y f + sin 6c R sin 9y,\ 
Set equations 44 and 45 equal and solve for Cy, leading to 

ar t Rmag (sin 9 a - tan 6c R cos 0$) 


Cy = 


(CTV\ tan 9 Cr + CTV 2 ) 


where 


CTV\ 

CTV 2 

CJT\ 

CT 2 

ct 3 

ct 4 

ct 5 

ct 6 


(CTi cos 8 Cv + CT 2 sin 8 Cv + CT 3 ) 

(' CT a cos 9 Cv + CT 5 sin 9 Cv + CT 6 ) 

VmagVf COS (j>y f sin 8y f sin 9y + R mag g cos (j)g COS 9g COS 9ft 
- VmagVf cos (j)y f sin 9y f cos 9y + R mag g COS (j) g sin 9g cos 9^ 

&TfRmag COS 9 p> 

VmagVf cos (f>y f cos 9y f sin Oy - R mag g cos 4> g cos 9 g sin 9p, 

— VmagVf cos <j)y f cos 9y, cos 9y — R m agg cos (f> g sin 9 g sin 9p, 
ClT f Rmag Sin 9 ff 


Now substitute equation 46 into equation 44 to find Cr , resulting in 

a Tf v mag (sin 9 y cos 9 Cv - cos 9y sin 9 Cv ) 


C R = ~- 
Define the following variables: 

CTV 3 

ct 7 

ct 8 

ct 9 

CT W 

Equations 46 and 48 can be written as: 

Cr = 

Cy = 


cos 8 Cr {CTV I tan 9 Cr + CTV 2 ) 


= CT 7 cos 8c v + CT 8 sin 9c v 
(It f Vmag sin 9y 
®Tf Vmag COS 9y 

— ClT f Rmag COS 9 ff 

— (XT f Rmag Sin 9 ff 


ctv 3 

cos 9q r (CTV i tan 9c R 
CT 9 tan 9c r + CT W 


CTV 2 ) 


( 45 ) 


(46) 


(47) 


(48) 


(49) 


(50) 


CTVi tan 9 Cr + CTV 2 

Comparison of Pi from equations 33, 34, 35, and 36 shows that a general equation for P can be written as: 


P = 


( Cr ■ V) 


— dTf 


D — 1 1 — 


D + 2 Et q + Ft 2 0 - 1 


(51) 


In Thrust Plane notation, this equation becomes 

CnYmag (cos 9 Cr cos 9y + sin 8 Cr sin 9y) 
/ 

= (lT f (Cy — 1 ) — dT 0 


(ylcfy + 2 C r C v <Vt 0 + C%r* - l) 


( 52 ) 
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Substitute equation 31 for $ into equation 52, isolate the square root term, square each side, and collect 
terms multiplied by C R and Cy. After some manipulation, we obtain 


0 = 


V^an (cos dp + sin dp tan d c J 2 - 


cos 2 0 Cb 


cos 2 0 Cr C 2 r 


+ (a Rf - a| o ) Cy 
- [ 2 a Tf Vmag (cos dp + sin dp tan 9 Cr 


+2 a\ o (cos 0c v + sin 9 Cv tan 0 Cr ) t 0 \ cos 9c r CrC v 
+ 2 (a Tf ~ a To ) Knag (cos dp + sin dp tan 0 Cr ) cos 9 Cr C r 
— ‘Zo-Tj ( ClTf ~ a To) Cy 

+ (ar f — a Ro ) (53) 


Note the following trigonometric identity: 42 


1 

cos 2 9 Cr 


(1 +tan 2 9 Cr ) 


(54) 


After substituting this identity and some manipulation, equation 53 can be written as 

0 = (CAi tan 2 0c R + CA 2 tan 9c R + CA 3 ) cos 2 9c r C r 
+ ( CTV 4 tan 9 Cr + CTV 5 ) cos 9 Cr C r C v 
+ ( CCi tan 9c R + CC 2 ) cos 9c r C r 

+ C D ± Cy + CD 2 C v + CD 3 (55) 

where 


CTV 4 

ctv 5 

CA 3 

CA 2 

ca 3 

CBj. 

CB 2 

cb 3 

CCi 

cc 2 

CD\ 

cd 2 

cd 3 


(CBi sin dev + CB 2 ) 

(CBi cos 9 Cv +CB 3 ) 

Vmag S^ 2 dp - 0%-J 2 

ZV^ag sin dp cos dp 

Vmag COS @V ~ a T Q T o 
— 2 ci r j' o i~o 
2 Vmag Sin Oy 
2 ClTf Vmag COS Oy 

2 {cL Tf &T 0 ) Vmag Sin Oy 

2 (aTf — ar o ) Vmag COS Oy 

( a Tf ~ a T 0 ^j 

—2aT f ( ar f — ar 0 ) 

(' a Tf -a To ) 2 


(56) 


The next step is to substitute C R and Cy from equation 50 into equation 55. After some manipulation, 
a quadratic equation for tan 9c R can be written as 

0 = Ctcr 2 tan 2 9c R + Ctcri tan 9c R + Ctcro (57) 
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where 


Ctcr'i = CTV*CA l + CTV?CD 3 

+ CTV [ CTV :i CC l + CTV3CTV4CT9 
+ CTV\C D 2 CTq + CD { CTg 
Ctcri = CTV3CA2 + CTV 3 (CW4CT-10 + CTV 5 CT 9 ) 

+ CTV t ( CD 2 CT w + CTV 2 CD 3 + CTV 3 CC 2 ) 

+ CTV-2 (CD 2 CT 9 + CTV\ cd 3 + CTV 3 CC\) 

+ 2CT>iCT 9 CTio 

Ctcro = CTV 3 CA 3 + CTV 2 C D 3 

+ ctv 2 ctv 3 cc 2 + CTV 3 CTV 5 CT w 

+ CTV 2 CD 2 CT 10 + CDiCTi 0 (58) 


Note that the coefficients in equations 58 are functions of 9c v and r D . 

A numerical search can be performed to solve for the values of 9q v and r G . For a given value of 9c v 
and r 0 , first assume a thrust sequence in order to set the initial and final values of the thrust acceleration 
magnitude. Now, compute all the values of CA i , CBi , CC), and CD, from equations 56. Also, compute the 
values of the CTi coefficients from equations 47 and 49. Note that all of these coefficients are only functions 
of t 0 . Now compute the values of the CTVi coefficients defined in equations 47, 49, and 56. The coefficients 
Ctcri-, Ctcri , and Ctcro defined in equation 58 can now be calculated. Two possible values for tan ^Cr 
can be computed from equation 57. From equation 54, two possible values of the absolute value of cos 9c R 
are given by 


I cos 0c R \ 


1 

y/l + tan 2 9 Cr 


(59) 


Using the two possible values for |cos0c< R | and tanf^R- two possible values for Cr and C v can be found 
from equations 50. The proper sign of cos 9c R is selected using the fact that Cr must be positive. With the 
values of tan 9c R and cos 9c R , two possible values of 9c R can be found. There is now enough information 
to calculate two possible values for Cr and Cy. The values of the integration constants P , G, and K 
can be computed from the equations defined in section VII. The V and R vectors can be computed from 
equations 27. Given the initial position, initial velocity, and gravity, the values of the final position and 
velocity can be computed from equation 26. The problem then reduces to finding the values of 9c v and r Q 
that meet the desired final position and velocity. The variable 9c v is bounded because it is an angle. It will 
be shown in section XI that the value of r Q can be bounded analytically. Thus, the problem is reduced to a 
two-dimensional bounded search. 


X. Reduction of One-Phase Optimal Problem to One Variable 


When the problem is constrained to a constant thrust acceleration magnitude, it is possible to analytically 
reduce the optimal problem to one variable. As in the previous section, it is assumed that a constant thrust 
direction solution is not valid. This is equivalent to the assumption that all of the vectors Cr,Cv, R, and V 
are non-zero and that none of them are parallel. For a constant thrust acceleration problem, the definitions 
of Pi, Gi, and K\ are given in equations 33. Set Pi = P, G 1 = G, K\ = K in equations 39 through 41, 
respectively. Define the ratio of Cr to Cy as 

P= % (60) 

Rearrange equations 39 through 43 in terms of p. Recall that the angle 9y is equal to zero by definition. 


CLt 


1 - V 1 + 2 P® t o + P 2 T% = pVmag COS 9c. 


(61) 


a-r In 


\A + 2 p$r 0 + p 2 tI + pr 0 + $ 


1 + 4 > 


pVmag ®hi 9q r 

(cos 9 C R sin 9 Cv - sin 9 Cr cos 9 Cv ) 


(62) 
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: a T r 0 \Jl + 2p^r 0 + p 2 T% 


— pRmag (COS 0 (J R COS 9 j_, si 1 1 0q r sill 0 p ) Vmag COS 0C\ 


0 = pRmag (cos 9 Cr sin 9ft - sin 9 Cr cos 9^) - V mag sin 9 Cv 
1 


(63) 

(64) 


0 = CiTf 


C v 


- 1 


+ pVf cos <t>v f (cos 9q r cos 9y, + sin 9q r sin 9 Vf ) 
+ g cos (j> g (cos 9c v cos 9 g + sin 9c v sin 9 g ) 


(65) 


Equations 61 through 64 represent a system of four equations in four unknowns: p , 9c R , 9 q v1 and r 0 . 
The variable Cy can be freely chosen to meet the condition in equation 65, but its value is not needed to 
solve the system of equations. Fill noted this fact: for a constant thrust or constant acceleration problem, 
the dimension can be reduced by ’’scaling the primer” vector 12 (p. 19). Multiply equation 61 by |t 0 and 
rearrange to obtain 

xOTToa/I + 2 p$T 0 + p 2 r 2 = -T 0 


2 x ■ u v - . -/ — ■ u ' r ■ o 2 U P^mag COS 

Setting the right-hand side of equation 66 equal to the right-hand side of equation 63 yields 

1 


(66) 


dT — pVmag COS 9c B 


1. 


= - pRmag (cos 9 Cr cos 9 R + sin 9 Cr Sin 0 k ) - -V mag cos 9 Cv 
Solving equation 67 for p and writing it in terms of tan 9c R results in 


P = 


Vmag COS Ocy 4“ ^T^~o 


cos 9 Cr 


^Rrriag sill 9 pj tail 9c R T ^V m ag^~o 2 R 


L mag 


cost 


Also solve equation 64 for p and write it in terms of tan 9c R to obtain 


V r 


P = 


mag sin 9c v 


cos 9c R R ma g (sin 9p — tan 0c R cos 9p) 

Setting equation 68 equal to equation 69 and rearranging yields 

0 = - [D Cr sin 9 Cv + E Cr cos 9 Cv + F Cr ] tan 9 Cr cos 9 Cr 
+ [ A c R sin 9 Cv + B Cr cos 9 Cv + C Cr ] cos 9 Cr 

where 


(67) 


(68) 


(69) 


(70) 


a c r 

— ^RmagVmag COS 9 ^ 

b c r 

— RmagVmag sin 9 

Cc R 

= arRmagTo sin 

d c r 

— 2 RmagVmag Sin 9 

e c r 

— RmagVmag COS 9 ^ 

f c r 

= CLt R mag^o COS 9 


(71) 

One possible solution is cos 9c R = 0. However, from equation 61 this would mean that SS = + 2p$r 0 + p 2 r 2 

1. This is equivalent to the special case when an off boundary thrust acceleration can solve the problem. 
This leads to a constant thrust direction solution. It was assumed at the beginning of the section that this 
case is not valid for the current derivation. Thus, for cos 9c R / 0 we have 

A Cr sin 9 c v + B Cr cos 9 Cv + C Cr 


tan 9 Cr = 


Dc r sin 9c v + Ec r cos 9c v + Fc R 


(72) 
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Now return to equation 61. Isolate the square root and square both sides to obtain 


0 = ^2 a 2 T <5>T 0 + 2 a T V mag cos 6 Cr ) P + (a| T o - ag cos 2 9 Cr ^ P 2 (73) 

One possible solution to this equation is p = 0. This implies that Cr = 0 and leads to a constant unit thrust 
direction with an off-boundary thrust acceleration magnitude. This solution is equivalent to dropping the 
constraint on final position. Thus, it is a special case that can be checked. Thus, for p ^ 0 we have 

_ _ 2 q|<f>T 0 + 2a T Vmag COS 9 Cr , , 

4b 2 - y^ag COS 2 9 Cr 

Set equation 74 equal to equation 69, substitute for the definition of $ from equation 31, use the trigonometric 
identity in equation 54, and rearrange to find a quadratic equation in tan 9q r ■ 

0 = [Cso sin 0c v + Cco cos 0 Cv + Co] 

+ [C S i sin 9 Cv + C C i cos 9 Cv + Ci] tan 9 Cr 

+ [C S 2 sin 9 Cv ] tan 2 9 Cr (75) 


where 


Cso 

2 _2t> t/3 

— LLrpT 0 Vmag v mag 

Cco 

— ‘2ci'jiT 0 R rna g sin. Or 

Co 

= 2cL'jpR rna g Vmag sin Or 

Csi 

— 2a T T 0 R rnag sin Or 

Cci 

= 2Q J rpT 0 R rna g COS Or 

Or 

= ^ClTRmagVmag COS 0 r 

C S 2 

= Vmag 


Now, substitute tan 9q r from equation 72 into equation 75 to obtain 

Ds sea sin 3 9c v 
+ Dsocs cos 3 9c v 
+ D S 2 C i sin 2 9 Cv cos 9 Cv 
+ D s i C 2 sin 9 Cv cos 2 9 Cv 
+ D S 2 C o sin 2 9c v 
+ -Ds°c 2 cos 2 9q v 
+ D s i C i sin 9 Cv cos 9 Cv 
+ D s i c osin9 Cv 
+ D S oci cos 9c v 

+ Dgoco = 0 


(76) 


(77) 
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where 


E>s 3 c° = CsoDc r + CsiAc r Dc r + CsiA? Cr 
Ds°c 3 = CcoEc r + CciBc r Ec r 
Ds^c 1 = CcoDc r + 2CsqDc r Eg r + CciAc r Dc r 
+ Csi {Ac r Ec r + Bc r Dc r ) + 2 Cs2Ag r Bc r 
Bs 1 c 2 = CsoEc r + 2CcqDc r Ec r + Cs\Bg r Ec r 
+ C C i (A Cr E Cr + B Cr D C r ) + C S 2Bc r 
E>s 2 c° = C 0 Dc r + 2CsqDg r Fc r + C\Ac R Dc R 

+ C S i ( A Cr F Cr + Cc r D C r ) + 2C S 2Ac r C Cr 
D soc 2 = C 0 E 2 Cr + 2CcoEc r Fc r + C\Bc R Ec R 
+ Cci ( B Cr F Cr + c Cr e C r ) 

B > s 1 c 1 — 2 C 0 Dc r Ec r + 2CcoDc r Fc r + 2 CsoEc R Fc R 

+ C 1 ! ( A Cr E Cr + B Cr D C r ) + Cci ( a Cr f Cr + c Cr d C r ) 

+ Csi (B Cr F Cr + Cg r E C r ) + 2C S 2 B Cr Cc r 
D s^c° = CsqEc r + 2C 0 Dc r Fc r + C\ ( Ac r Fc r + Cc R E>c R ) 

+ CsiCc r Fc r + Cs2Cc r 

Ds°c 1 = CcoFc r + 2C 0 Ec r Fc r + C\ ( Bc R Fc R + Cc R Ec R ) 

+ CciCc r Fc r 

D S 0 C 0 = C 0 F 2 r +CiCc r Fc r (78) 

It is shown in reference 24 that all coefficients that do not multiply with a power of sinf?cv are zero, so 
Dsoc 3 = Dgoc 2 = Dsoci = Dsoc° = 0. Equation 77 then reduces to 

D S 3c° sin 3 0c v 
+ D S 2 C i sin 2 0 Cv cos 9 c v 
+ D s i c 2 sin 9 c v cos 2 9 Cv 
+ Dg 2 G° sin" 9 Gy 
+ D S i C i sin 9 Cv cos 9 Cv 

+ D s i c osm9 Cv =0 (79) 

One possible solution is sin$cv = 0- This means that Cy would be aligned with V. From equation 69, it 
can be seen that p = 0 when sin^cy = 0, and when p = 0, then Cr = 0. As discussed previously, it has 
been assumed for purposes of this derivation that this is not the case. Thus, for sin0<~v 7^ 0, equation 79 
reduces further to 

D s 3 c° sin 2 9 c v 
+ D s 2 C i sin 9 c v cos 9 Cv 
+ D S ic 2 cos 2 9 c v 
+ Ds 2 c° sin 9 c v 
+ D S i c i cos 9 c v 
+ Ds^c° = 0 

Note the following well known trigonometric identity: 42 

sin 2 9c v + cos 2 9 c v = 1 
=> cos 2 9c v = 1 — sin 2 9c v 

=> cos 9 g v = ±\J (l — sin 2 9c v ) 


(80) 


(81) 
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Use this identity in equation 80, rearrange to isolate the square root term, square both sides, and rearrange: 


(Ds sc 0 — D S i C 2 ) 2 + D 2 S 2 C i sin 4 9 Cv 


+ [2 (D S 3c° - Ds^c 2 ) D S 2 c° + 2 D S 2 C iD s i c i] sin 3 9c v 

+ [2 ( D S 3c° — D S ic 2 ) (D S ic 2 + Ds^c 0 ) + D 2 S 2 C o + D 2 S i c i — D 2 S 2c i\ sin 2 9c v 
+ [2 ( D S ic 2 + D S i C o) D S 2 C o — 2D S 2 C iD S i C i] sin 9c v 


( D S ic 2 + D S ic °) 2 — D 2 s i c i 


= 0 


(82) 


Note that the coefficients of equation 82 are functions of only r 0 . 

A numerical search can be performed to solve for the value of r 0 . For a given value of r 0 , compute all 
the values of Ac R , Bc R , Cc R , Dc R , Ec R , and Fc R from equation 71. Compute the values of Cgo , Cco , Co, 
Csi, Cci, Ci, and Cg i from equation 76. Compute the values of Dg3Qo, Dg 2 Qi, Dg \Q 2 , Dg 2 qo, DgiQi, and 
DgiQO from equation 78. The roots of the quartic defined in equation 82 can now be computed to give four 
possible values of sin^cy- Using the inverse sine function, it is possible to find four possible values of 9c v 
between — 7t/2 and n/2. However, it is possible for the proper value of 9c v to fall between — 7r and 7r. Thus, 
from the four possible values of sin(?cy, there are eight possible values of 9c v - The proper values of 9c v 
can be resolved by checking each of the eight possible values in equation 80. This will result in four possible 
values of 9c v between —i r and n. Equation 72 can be used to compute the possible values of tan 6^- As in 
the previous section, equation 54 can be used to find the absolute value of cos 9q r : 


| cos 9 Cr | 


1 

\/l + tan 2 9 Cr 


(83) 


The possible values of p are calculated from equation 69. The proper sign of cos 9q r is selected using the 
fact that p must be positive. With the values of tan 9c R and cos 9c R , the values of 9c R can be found. The 
proper value of Cy can be computed from equation 65. There is now enough information to calculate four 
possible values for Cr and Cy. The values of the integration constants P, G, and K can be computed from 
equations 33. The V and R vectors can be computed from equations 27. Given the initial position, initial 
velocity, and gravity, the values of the final position and velocity can be computed from equation 26. The 
problem then reduces to finding the values of r c that meet the desired final position and velocity. It will be 
shown in section XI that the value of r Q can be bounded analytically. Thus, the problem is reduced to a 
one-dimensional bounded search. 


XI. Bounding the Optimal Control Problem 

XI. A. Optimal Control Problem Lower Boundary 

A lower boundary on the performance index of the problem can be found by solving the problem with no final 
position constraint. In this case, the optimal control is a constant thrust direction with a thrust acceleration 
which can be off its lower and upper boundaries. This subproblem can be solved analytically. Keeping the 
final velocity constraint and dropping the final position constraint from the problem can be achieved by 
choosing ctr = 0. This results in the following condition: 

Cr = 0 (Inactive final position constraint) (84) 

Ay = Cy = constant (Inactive final position constraint) (85) 

From equations 14, it can be seen that E = F 
reduces to 

u = 

The unit thrust direction is constant for the entire optimal trajectory. The thrust acceleration magnitude 
switching structure (SS) was defined in equation 18. Since both Ay and u are constant, then the SS is 
constant over the entire optimal trajectory. Thus, the thrust acceleration magnitude will be either on its 


= 0. From equation 15, the optimal unit thrust direction 
Ay 


V 


Gy 

\/D 


(86) 
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upper boundary, lower boundary, or off-boundary for the entire optimal trajectory. Since the unit thrust 
vector is constant, equation 26 can be written as: 




ot dr 


u 




clt dr dr 


u 


(87) 


The integrated equation for velocity can be used to determine the constant unit thrust acceleration direction 
required to meet the final velocity constraint as 


(V f - V 0 ) - gr 0 

fo° a T d A 


(88) 


This shows that the optimal unit thrust direction is parallel to V. The unit thrust vector must have length 
of one. Taking the vector dot product of equation 88 with itself and setting it equal to one results in 



[(<?-s)K 2 - 2 (v f -Vo) -9 


To + 


(Vf-v 0 )-(v f -v 0 ) 


(89) 


Equation 89 gives the value of the performance index, the integral of the thrust acceleration, as a function 
of the time-of-flight. The average value of the thrust acceleration is given as 


CLj 1 — 


/;• ardr _ ^-9)\ - 2 [{?, ~ <t) • <?] + [(?/- ft) - (v, - V.) 


(90) 


Since the thrust acceleration can vary between aT min and OT max , then the average value of thrust acceleration 
over any time interval is also bounded by the maximum and minimum acceleration values. 

Consider the case when the thrust acceleration is off-boundary. An equation for the time-of-flight that 
minimizes the unbounded performance index can be written by taking the derivative of equation 89 with 
respect to r 0 and setting it equal to zero. This is referred to as the “unbounded performance index” because 
it does not account for the thrust acceleration bounds, hence 


(Vf-V 0 )-g 

(9-9) 


(91) 


It is simple to verify that the unbounded time-of-flight will cause the unit thrust direction to be perpendicular 
to the gravitational vector, as shown in equation 22. It will also make SS=1. Substitution of equation 91 
into equation 89 results in the unbounded optimal performance index as 



axdr 


J unbounded 



(Vf-Vo) -9 

(<?•(?) 


(92) 


The average thrust acceleration necessary to achieve this performance index value can be found by substi- 
tution of equation 91 into equation 90 and use of the trigonometric identity in equation 54. 


= ls|tan<W 9 (93) 

where 6/yv g is the angle between the gravity vector, g , and the vector representing the change in velocity 
from the initial to the final state, (Vf — Vo'j. Note that this angle is bounded such that 0° < 9a V g < 180°. 

Recall that in this special case, the thrust acceleration is applied perpendicular to the gravity vector. 
When tan OAVg < 0, this implies that the problem is not controllable when thrust is applied only perpendic- 
ular to the gravity vector. When this occurs, the thrust acceleration magnitude must be on the maximum 
or minimum boundary. 
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For the case when 90° < 0AVg < 180°, the change in velocity requires the thrust to act against gravity. 
Thus, a component of the acceleration must be in the direction parallel to gravity. In order to minimize 
the gravity losses, the minimum time solution is sought. This means that the optimal solution uses the 
maximum acceleration possible. 

For the case when 0° < 6a V g < 90°, the change in velocity allows some use of gravity to achieve it. The 
optimal component of acceleration in the channel parallel to gravity is zero. In this case, gravity is allowed 
to carry the vehicle along the gravity channel while all control effort is focused on the velocity perpendicular 
to the gravity vector. However, it is possible that this “unbounded” optimal solution results in an average 
acceleration which is outside the acceleration limits. In this case, the acceleration limit closest to the average 
acceleration is the bounded optimal solution. From equation ??, a quadratic equation for the time-of-flight 
when the thrust acceleration is on one of its boundaries can be written 


[(9- 9) 4 lim J r ° 2 - 2 [(V f V 0 ) ■ g] t 0 + [(V f - V 0 ) ■ (v f ~ V 0 ) 


= 0 


Only positive roots of this equation are possible solutions. The value of the performance index is then given 
by the value of multiplied by the positive roots of this equation. For real world problems, only real 

solutions of this equation are allowed. By examining the discriminant of this quadratic equation, a necessary 
condition on a,T limit for real solutions is found to be: 24 

a T H mu - I <7 1 2 sin2 (94) 


Thus the physics of the problem impose a lower limit on the possible thrust acceleration magnitude. For a 
solution to this problem, the maximum thrust acceleration limit must be greater than this value. 

There is now enough information to find the solution to this subproblem. The equations derived in this 
section can be used to evaluate the performance index for the cases when a-r = aT min , cit = «T r „ ox . , and 
a T mirl < clt < ciT max ■ The solution which gives the minimum value of the performance index while satisfying 
the constraints is a lower bound of the problem. 


XI. B. Optimal Control Problem Upper Boundary 


Any feasible solution to the problem is an upper bound on the performance index. Consider a trajectory 
that is comprised of two segments with constant thrust acceleration magnitude and direction. It is desired 
to find the thrust magnitude, thrust direction, and segment times that will take a vehicle from an initial 
position and velocity state to a final position and velocity state. 

The equations of motion from equation 1 can be easily integrated along each segment, due to the constant 
thrust acceleration and gravity. The first segment runs from the initial state and time to the unknown 
switching state and time. The second segment runs from the unknown switching state and time to the final 
state and time. Combining the two equations to eliminate the unknown switching state and time results in 


where 


V 

R 


[a T2 T 2 ] u 2 + [arUi] ui 


, a T 2 T 2 


U 2 + 


-a Tl (r 2 + r) 7~i 


Ui 


(95) 


T 

Tl 

T2 


(tf — t 0 ) = Total time-of-flight 

(t s — t 0 ) = Time-of-flight for segment 1 

(tf — t s ) = Time-of-flight for segment 2 


(96) 


Note that 


r = n + r 2 


(97) 


Define the ratio of t 2 to r as 


T2 

P = — 


(98) 


T 

Note that since 0 < t 2 < r, then it must be true that 0 < p < 1. Now t\ and t 2 can be written in terms of 
r and p as 


n = r(l -p) 

r 2 = rp (99) 


21 of 32 


American Institute of Aeronautics and Astronautics 



Substitute equations 99 into equations 95 to obtain 


V 

— = [CIT 2 P\U2 + [UT, (1 — P)]U1 

T 

R 

^2 — 


ar 2 P 


u 2 


(1 - P 2 ) 


Ui 


Solving these equations for the thrust direction unit vectors yields 


(100) 


u 1 


U2 


2 R - Vrp 
a Tl r 2 (1 - p) 

- 2 R + Vt (1 + p) 

ClT 2 T 2 p 


( 101 ) 


Note that the thrust direction vectors must be unit vectors by definition. Taking the vector dot product of 
Mi and u 2 from equation 101 with themselves, setting them equal to one, and rearranging results in 


a 2 Ti 


a T 2 


4 (R ■ Rj - 4 (R ■ V^j t P +(v ■ V'j t 2 p 2 
r 4 (1 - p) 2 

4 [r ■ 1?) - 4 (r ■ V^j T (1 + p) + (v ■ V-) r 2 (1 + pf 


(102) 


For any given value of r, it will be useful to know the value of p which makes the thrust acceleration of each 
segment equal to a given value (for example, the maximum or minimum acceleration limits). Rearranging 
equations 102 results in quadratic equations for p as a function of ar and r. For ctr, , we have 


For ar 2 , we have 


Two quadratic equations 


+ 

+ 


-a 2 Tl r 4 + (V ■ V^j t 2 
24 1 t 4 -4(R-F)t 
-a 2 T j A +4 (r- Rj 


I P 

= 0 


2 4 

-a T2 r 


(v-v^t 2 p 2 
+ 2 (v ■ y) r 2 - 4 [r ■ v^j t p 

+ (v ■ V^j T 2 - 4 (r • v'j T + 4 (r ■ R^j 


= 0 


a 2 x 2 + a±x + a 0 = 0 

b 2 x 2 + b\X + b 0 = 0 


(103) 


(104) 


(105) 


can be solved by the same value of x only when the following condition holds: 24 

0 = [a 2 5i - CI 1 & 2 ] [a\b 0 - a 0 b{\ - [a 2 b 0 - a 0 b 2 f 


(106) 


This condition can be used with equations 103 and 104 to find a polynomial equation for r as a function 
of the initial state, final state, gravity, and thrust acceleration magnitude on each segment. The following 
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condition is obtained: 


0 = { [(V- F)r 2 -a^r 4 -4 (r ■ V - ) r + 2 (v • V") r 2 ... 

— — 4 ^7?. • r + 2o^ 1 r 4 ^y • r 2 — a^r 4 |... 

* | -4 ■ v) t + 2oy x r 4 (y • y) r 2 - 4 (r ■ V^j T + 4 (r ■ r) J ... 

- a(r-R) -a 2 Ti r 4 -4(i?-y)r + 2(y-y)r 2 }... 

- { (v • y) r 2 - a^r 4 (y.y)r 2 -4^.y)r + 4(^.^)]... 

- a(r-R) -a 2 Ti T 4 (v ■ V^j T 2 - a 2 T 2 T 4 | 2 (107) 

Each coefficient of the polynomial is a function only of the initial state, final state, gravity, and the thrust 
acceleration magnitude on each segment. An initial attempt to collect terms of r by hand was abandoned. 
Instead, a Matlab script was written to symbolically convolve each subpolynomial in tau in order to define 
the coefficients. The resultant 10th order polynomial in r is much too large to put into print. 

There is now enough information to solve the problem for a given initial state, final state, gravity, and 
maximum and minimum thrust acceleration. The coefficients of equation 107 are evaluated, assuming thrust 
acceleration magnitudes for each segment. All positive, real roots of this equation are potential solutions for 
r. If no positive, real roots exist, different thrust acceleration magnitudes are chosen. Once potential values 
of t are available, the roots of either equations 103 or 104 are found. The roots that are between zero and 
one are potential values of p. The value of the performance index is given by: 

a T dr = a^Ti + clt 2 t 2 = [a^ (1 — p) + clt 2 p\ t (108) 

If multiple solutions exists, the one with the lowest performance index is chosen. The unit thrust direction 
for each segment can then be computed from equations 101. 

XII. Numerical Solution of Reduced Problem 

An algorithm to solve the optimal control problem with a constant thrust acceleration magnitude was 
developed. This algorithm is a mechanization of the equations from section X. Lower and upper bounds 
on the time-of-flight are given in section XI. The algorithm is a one-dimensional bounded search developed 
specifically for the problem at hand. 

In order to test the proposed guidance algorithm in various scenarios, a simulation was developed in 
Matlab. 43 The simulation models a vehicle as a point mass over a rotating, spherical planet with an inverse 
gravitational field. Both the navigation and control systems are assumed to be perfect. The guidance is 
called at a specified rate, and the solution is assumed to be instantaneously available. The vehicle rocket 
engine is assumed to be perfectly expanded and throttleable. In practice, the guidance commands a thrust 
acceleration level. Both nominal runs and monte carlo runs are completed. Each monte carlo run includes 
100 dispersed trajectories. For details on the models of the simulation, see reference 24. 

A numerical optimizer was used to test the optimality of the solution from the simulation. The numerical 
optimization uses the Legendre pseudospectral collocation method. A pseudospectral differentiation matrix 
is used to write the equations of motion as a set of nonlinear algebraic equations. The problem is then 
written numerically as a parameter optimization problem, where the states, controls, and time-of-flight are 
manipulated by a numerical optimizer to achieve an optimal solution. Details of the optimization method 
can be found in reference 44. 

XII. A. Test Case Definition 

The lander begins in a circular orbit with an altitude of 100 kilometers and an inclination of 60 degrees. It 
is at its maximum southern latitude when a deorbit burn is completed to put the lander into an elliptical 
transfer orbit with an apolune altitude of 100 km and a perilune altitude of 20 km. When the lander 
reaches the perilune of the transfer orbit, the descent phase begins with a braking burn aligned opposite 
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the planet relative velocity vector. This portion of flight is not modeled in the simulation. When the 
altitude drops below a threshold of 6 km, the closed-loop guidance commands the thrust acceleration vector 
to achieve a desired final position and velocity state. Guidance is called every 10 seconds to generate new 
thrust commands. Between guidance calls, the thrust acceleration command profile from the last guidance 
call is followed open loop. In order to simulate a possible hazard avoidance maneuver, the target may be 
redesignated at a total range-to-target of 2 kilometers. If this occurs, the new landing target is loaded into 
guidance at the time of redesignation; the guidance has no apriori knowledge of the new landing site. When 
the time-to-go of the optimal flight phase drops below 10 seconds, guidance is no longer called, and the final 
guidance commands are followed open loop. At the termination of the optimal flight, the vehicle is over the 
landing site with a very small negative altitude rate and no horizontal velocity. The pitch up and vertical 
landing phases are not modeled in the simulation. 

The guidance algorithm computes commands in a topocentric guidance frame. For this frame, the X-axis 
is East, the Y-axis is North, and the Z-axis is Up. The guidance algorithm assumes that all gravitational 
acceleration is directed along the negative Z-guidance axis. During each call to guidance, a corrected gravity 
term is added to account for the centripetal acceleration of the topocentric frame. 


0 


0 

0 

+ 

0 



V? 

1 

O 


* horizontal 

Rmag 


(109) 


where 


g 0 = Gravitational acceleration magnitude at surface of planet 
Vi horizontal ~ Inertial velocity in the local horizontal plane 
Rmag = Radial position magnitude 


The initial state of the lander is given in Table 1. The selected landing targets are given in Table 2. The 
lander vehicle properties are given in Table 3. The properties of the moon are given in Table 4. The monte 
carlo dispersions are given in Table 5. 


Variable 

Units 

Initial State 

Range North of Target 

m 

3100 

Range East of Target 

m 

-21120 

Geocentric Altitude 

m 

6000 

Planet Relative Velocity Magnitude 

m/s 

350 

Planet Relative Topocentric 
Flight Path Angle 

deg 

-18.65 

Planet Relative Topocentric 
Azimuth Angle 

deg 

104 


Table 1. Lunar Test Case Initial State 
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Variable 

Units 

Landing Site 

Primary 

Divert 

Geocentric Latitude 

deg 

58.9 

58.91649 

Longitude 

deg 

146.73 

146.76298 

Geocentric Altitude 

m 

100 

100 

Geocentric Altitude Rate 

m/s 

-5 

-5 

Planet Relative 
Horizontal Velocity 

m/s 

0 

0 


Note: The divert landing site is 0.5 km north and 1.0 km east 
of the primary landing site. 


Table 2. Lunar Test Case Landing Targets 


Variable 

Units 

Value 

Thrust Acceleration Magnitude 

m/s 2 

5.5 

Gravitational Acceleration Magnitude 
at Surface of Planet 

m/s 2 

1.635 


Table 3. Lunar Test Case Vehicle Properties 


Variable 

Units 

Value 

Equatorial Radius 

km 

1737.4 

Gravitational Constant 

km 3 /s 2 

4.9028xl0 3 

Sidereal Rotation Period 

hours 

655.728 


Table 4. Planetary Properties of the Moon 


Variable 

Units 

Mean 

Uniform 

Dispersion 

Range North of Target 

m 

3100 

±500 

Range East of Target 

m 

-21120 

±500 

Geocentric Altitude 

m 

6000 

±100 

Planet Relative 
Velocity Magnitude 

m/s 

350 

±5 

Planet Relative Topocentric 
Flight Path Angle 

deg 

-18.65 

±0.25 

Planet Relative Topocentric 
Azimuth Angle 

deg 

104.0 

±0.25 


Table 5. Lunar Test Case Monte Carlo Dispersions 
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XII. B. Flight to the Primary Target 

Table 6 contains results for flight to the primary target. A comparison is given between the optimized 
trajectory, the nominal trajectory from the simulation, and statistical results of the 100 dispersed trajectories 
at the final flight time. Target accuracy for all cases is very good, and the performance index from the nominal 
simulation is very close to the numerically optimized performance index. The performance index shows a 
variation of roughly 90 meters/second due to the dispersed initial state in the monte carlo analysis. 


Parameter 

Units 

Target 

Optimizer 

Nominal 

Min 

Mean 

Max 

Total Range 

m 

0 

1.25e-l 

8.55e-2 

8.40e-5 

8.99e-2 

2. 01e-l 

North Range 

m 

0 

1.78e-2 

-8.55e-2 

-1.99e-l 

-8.48e-2 

7.53e-3 

East Range 

m 

0 

-1.24e-l 

0.00 

-6.34e-2 

1.12e-2 

7.32e-2 

Geocentric 

Altitude 

m 

100 

100.006 

101.217 

100.538 

101.393 

102.060 

Geocentric 
Altitude Rate 

m/s 

-5 

-5.000 

-4.844 

-4.899 

-4.833 

-4.782 

Planet Relative 
Horizontal Velocity 

m/s 

0 

1.465e-3 

1.26e-2 

3.80e-4 

1.123e-2 

1.841e-2 

Time-of-Fliglit 

s 

- 

75.2567 

75.3535 

74.706 

78.896 

91.257 

Performance Index 

m/s 

- 

413.9118 

414.4443 

410.883 

433.928 

501.914 


Table 6. Primary Target: Comparison of Optimizer and Simulation at Final Time 


Figures 1(a) through 1(h) show trajectory plots for flight to the primary target. For each plot, the solid 
red line shows the nominal trajectory from the simulation, the dashed cyan line shows the trajectory from 
the optimizer, and each solid blue line represents one dispersed trajectory. From these figures, it can be seen 
that the nominal trajectory matches the optimized trajectory very well. It can also be seen that each monte 
carlo flight hits the desired final target state. 

Figure 1(a) shows the geocentric altitude versus the total range to the target. Figure 1(b) shows the 
time history of the planet relative velocity magnitude. Of particular interest is Figure 1(c), which shows the 
groundtrack in terms of geocentric latitude versus longitude. Each dispersed trajectory indeed hits the final 
target point, but 65% of the cases overfly the target by more than 10 meters before turning around. This 
can also be seen in Figure 1(d), which show the time history of the planet relative topocentric flight path 
angle. The variations seen in flight path angle near the end of the trajectories show the flight characteristics 
of a vehicle correcting for overflight. It can also be seen that the flight path angle varies slightly between 
the nominal and optimized trajectories near the end of flight. 

Figures 1(e) and 1(f) show the thrust direction angles. Differences can be seen between the nominal 
and optimized trajectories which lead to a slightly better solution for the optimizer. For the monte carlo 
results, the variations at the end of the trajectory show vehicle response to the overflight. However, it has 
been verified with the numerical optimizer that the overflight trajectories are indeed optimal for flight from 
each particular initial dispersed state. The overflight is predicted on each call to guidance, beginning with 
the very first call. This is due to the fact that a constant thrust acceleration magnitude, ax, is used for the 
entire trajectory. For a particular dispersed initial state, a different value for ax would minimize the amount 
of overflight. However, for the value of a-r used in this test case, the overflight is fuel optimal for many of the 
dispersed trajectories. No analysis was performed to optimize the thrust acceleration magnitude for each 
dispersed trajectory. 

Figure 1(g) shows the total angle of attack. The difference between the nominal and optimized trajectories 
is most pronounced in these plots. Note that the optimal total angle of attack is not zero, although it is 
small. Thus, a gravity turn which thrusts opposite the velocity vector is not optimal for the constraints 
of this test case. The same trend is seen in the monte carlo trajectories. Figure 1(h) shows the estimated 
time-to-go as a function of time for each trajectory. No line is given for the optimizer because it is flown 
off-line. As the guidance is called every 10 seconds, the time-to-go decrements nearly linearly with time for 
each trajectory. 
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XII. C. Flight to the Divert Target 


Table 6 contains results for flight to the divert target. For this scenario, the guidance is initially given the 
primary landing target. At a total range-to-target of two kilometers, the divert landing target defined in 
Table 2 is given to the guidance. The divert landing site is five hundred meters north and one kilometer 
east of the primary landing site. For direct comparison, the exact initial state dispersions used for the flight 
to the primary target are used here. A comparison is given between the nominal trajectory to the primary 
target, the nominal trajectory to the divert target, and statistical results of the 100 dispersed trajectories at 
the final flight time. Note that the target accuracy is very good. The performance index shows a variation 
of roughly 40 meters/second due to the dispersed initial state. 


Parameter 

Units 

Target 

Nominal Primary 

Nominal Divert 

Min 

Mean 

Max 

Total Range 

m 

0 

8.55e-2 

3.62e-l 

8.00e-3 

2.47e-l 

4.19e-l 

North Range 

m 

0 

-8.55e-2 

-2.58e-l 

-2.62e-l 

-1.71e-l 

8.00e-3 

East Range 

m 

0 

0.00 

2.55e-l 

0 

1.74e-l 

3.30e-l 

Geocentric 

Altitude 

m 

100 

101.217 

101.730 

100.540 

101.388 

101.854 

Geocentric 
Altitude Rate 

m/s 

-5 

-4.844 

-4.818 

-4.894 

-4.837 

-4.812 

Planet Relative 
Horizontal Velocity 

m/s 

0 

1.26e-2 

4.098e-2 

1.894e-3 

2.864e-2 

4.232e-2 

Time-of-Flight 

s 

- 

75.3535 

78.675 

77.543 

78.555 

84.546 

Performance Index 

m/s 

- 

414.4443 

432.7125 

426.487 

432.053 

465.003 


Table 7. Divert Target: Comparison of Divert and Primary Trajectories at Final Time 


Figures 1 (i) through 1 (p) show data for the Lunar Test Case nominal flight to the primary and divert 
landing sites. For each plot, the solid red line shows the nominal trajectory to the divert target, the dashed 
cyan line shows the nominal trajectory to the primary target, and each solid blue line represents one dispersed 
trajectory. The two nominal trajectories are exactly the same until the total range-to-target reaches a value 
of two kilometers. Figure 1 (i) shows the geocentric altitude versus the total range to the target. The 
discontinuity in the range occurs when the divert target is loaded into guidance. Figure 1 (j) shows the time 
history of the planet relative velocity magnitude. Figure 1(c) shows the groundtrack in terms of geocentric 
latitude versus longitude. Note that none of the trajectories overfly the divert target. Figure 1(1) shows the 
time history of the planet relative topocentric flight path angle. From these plots, it can be seen that each 
trajectory hits the desired final target state. 

Figures l(m) and l(n) show the thrust direction angles. The guidance response to the divert target 
can be seen clearly. Especially of interest is Figure 1 (o) , which shows the total angle of attack. Flight to 
the primary target has a small total angle of attack. This shows that the optimal thrust direction is very 
nearly opposite the velocity vector. However, for the divert trajectory, it can be seen that thrusting along 
the velocity vector is not desired. For both the nominal and monte carlo flights, it can be seen that the 
optimal thrust profile is not aligned opposite the velocity vector during the divert maneuver. The estimated 
time-to-go versus time is shown in Figure 1 (p) . A discontinuity occurs when the divert target is loaded into 
guidance. Note that except for this discontinuity, the estimated time-to-go is nearly linear with time. 

XIII. Conclusions 

Powered descent is a problem often characterized by a desire to safely achieve a specified landing target 
while minimizing the amount of fuel used in the maneuver. Integration of the equations of motion with 
the well known fuel optimal bilinear tangent steering law results in a system of five nonlinear equations in 
five unknowns. It was observed that the optimal unit thrust vector lies in a plane completely defined by 
the initial position, initial velocity, final position, final velocity, gravity, and the time-of- flight. The major 
contribution of this research is the use of this Thrust Plane to reduce the dimensionality of the nonlinear 
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problem. The problem with lower and upper thrust acceleration bounds can be reduced to a system of two 
equations in two unknowns: the time-of-flight and the angle of the velocity Lagrange multiplier in the Thrust 
Plane. The problem with a constant thrust acceleration can be reduced to one equation in one unknown: the 
time-of-flight. In addition, analytical lower and upper bounds on the time-of-flight have been derived. Thus, 
the problem is reduced to either a bounded two-dimensional search or a bounded one-dimensional search. 
The bounded one-dimensional search that defines the problem with constant thrust acceleration has been 
mechanized. It was demonstrated as a potential guidance algorithm for terminal descent to the surface of 
the moon. Its potential use for hazard avoidance divert maneuvers was shown. Several ares for future work 
have been identified: 

1. The bounded one-dimensional search for the constant thrust acceleration solution was mechanized 
efficiently. Work was begun on a similar mechanization of the bounded two-dimensional search for the 
problem with bounded thrust acceleration. However, an efficient algorithm has not yet been found. 

2. The reduction of dimensionality was done for problems with bounded thrust acceleration. However, 
some problems may be better described by a bounded thrust magnitude. The concept of the Thrust 
Plane may be applicable to such problems, allowing for reduction of dimensionality for bounded and 
constant thrust magnitude solutions. 

3. During the course of this research, it was noted that the bounded thrust acceleration problem could 
be reduced to one variable if enough numerical precision were available to compute the coefficients of a 
polynomial in time-of-flight with order greater than twenty. The coefficients themselves are high order 
polynomials of the parameters of the problem, and thus are prone to numerical error. This avenue was 
not investigated in the current research. 

4. The derivation of the optimal control law does not constrain the trajectory to fly above the surface of 
the planet. The trajectory is determined by the initial state, final state, and gravity. It is very possible 
that the optimal trajectory will fly subsurface. For the test cases of this research, the initial states 
were chosen to avoid subsurface flight. This is not difficult, but care must be taken to choose a proper 
initial state such that the trajectory will fly above the surface for all dispersed states. Future work 
may be done to include a constraint on altitude to ensure flight above the surface. 

5. For guidance implementation of the constant thrust acceleration algorithm, it may be desired to search 
for the best value of the constant thrust acceleration. This may allow for lower fuel usage for dis- 
persed trajectories. It may also be possible to select proper thrust acceleration magnitudes that avoid 
subsurface flight. 

6. The algorithm was tested in a three clegree-of-freedom simulation with perfect control and perfect 
navigation. While this is a good environment for proof-of-concept work, the algorithm should be 
tested in a high-fidelity, six degree-of- freedom environment. 
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